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A new improved transfer matrix method (TMM) is presented. It is shown that the method not 
only overcomes the numerical instability found in the original TMM, but also greatly improves the 
scalability of computation. The new improved TMM has no extra cost of computing time as the 
length of homogeneous scattering region becomes large. The comparison between the scattering 
matrix method(SMM) and our new TMM is given. It clearly shows that our new method is much 
faster than SMM. 
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Transfer matrix method (TMM) has been a useful ap- 
proach for studying the physical properties. There have 
been many publications on application of the TMM, such 
as studies of Ising model 0, 0, @, Q > quantum spin Q , 
electronic transport [f| 0] , and electronic state in quasi- 
periodic and aperiodic chains [1, |i| . TMM is also widely 
applied in studying the propagation of electric-magnetic 
wave [ioj . elastic wave [n| and light wave [l2[ in multi- 
layer systems. 

Original TMM(OTMM) can be efficiently used for the 
studies of periodic system, one boundary problem be- 
tween two homogeneous media and the scattering prob- 
lem of a small region sandwiched by two leads. How- 
ever, OTMM has a fundamental shortcoming due to its 
numerical instability when the number of transfer steps 
is beyond a critical value. This has been a major lim- 
itation to the application of OTMM. Some approaches 
have been developed to calculate the transport through 
a longer scattering area such as scattering matrix method 
(SMM) [l3|, [TJ, [l5| , extended transfer-matrix technique 
(ETMT)[tJ and Green's function approach(GFA)[l6j]. In 
SMM, the middle scattering region is cut into some 
smaller sections, then the scattering matrices of these 
sections are found by means of TMM and combined to- 
gether to get total scattering matrix(SM) recursively. In 
ETMT, they used a technique to cancel the increasing 
modes step by step for each slice to avoid the numer- 
ical instability. Although the numerical instability in 
SMM and ETMT does not exist anymore, CPU time 
cost grows linearly when the length of middle scatter- 
ing region increases. GFA too. Based on a continu- 
ous Schrodinger equation in electron wave guide without 
spin-orbit coupling (SO C)[l7| and with SOCQIj], the sta- 
ble solution of transport were also achieved. However, 
there are infinite number of evanescent modes that pose 
tremendous numerical difficulty, and the approximation 
of limit number of evanescent wave must be done. But in 
some case, the contribution of evanescent wave is signif- 
icant, even more in the presence of SOCfH, In this 
letter, we present a new improved TMM (NITMM) that 
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has overcome such non-physical numerical instability of 
OTMM. Our study will focus on the discrete version of 
Schrodinger equation with Rashba SOC[13| and will find 
the exact solution of electron transport with and without 
SOC. Our new method not only gives numerically stable 
solution at any length of the sample, but also provides 
high performance in computation. That is, no extra com- 
puting costs at increasing length of sample. Meanwhile 
the simplicity of OTMM is still maintained. 

General formulae of TMM. A 2d strip with the ge- 
ometry of a bar sandwiched by two semi-infinite leads is 
studied. The Hamiltonian of 2d electron gas is 

H = ^-^+V(x,y)+H so , (1) 
2m* 

where V{x, y) is the potential including confinement 
boundary, H so the remainder of the Hamiltonian that 
may contain spin-orbit coupling. After discretization of 
the Schrodinger equation H^(x, y) = E^(x, y), the strip 
is replaced by a lattice array with N infinite long chains. 
The lattice constant is assumed to be a. Wider strip needs 
larger chain number N to represent. The Schrodinger 
equation can be transformed into the following transfer 
matrix (TM) equations: 

nm 
Ti)$i, (2) 
i— l 

where $i = {4>i,N,t, 4>i-\,N,-[, 4>i,N,l> <t>i-i,N,l, ■ • • , 0i,i,T> 
<Pi-i, i, Tj 4>i, i)*- &i is a matrix with dimension 
AN x 1. The superscript t means transpose of matrix. 
The element 4>i,j,a m &i is t ne value of wave function 
ty(x,y) at site with spin index er. The lattice in- 

dices are | i £ (—00,00), j E (1,N)}. The left 

lead is in the region i £ (—00, — L — 1], the right lead in 
i G [£ + 1, 00), and middle bar in i € [— L + l, L — 1]. Two 
interfaces are at i = —L and i = L. % is a 47V x AN TM 
between the wave functions <&i and Two leads can 

be different or the same, and here we assume they are the 
same. Further, we assume that the leads and middle bar 
are homogeneous. Thus we have four different transfer 
matrices: T/ is the TM in two leads, T so in bar, T$ L and 
Ts R at left and right interface respectively. 

Fixing adjustable parameters like hopping constant to, 
on-site potential and the energy of electron E, we can ob- 
tain all elements of four TM matrices {Ti,T so ,Ts L ,Ts R }. 



2 



It can be verified that the determinants of Ti and T so 
within homogeneous regions are 1. Then two transform 
matrices Ui and U so can be numerically solved to diago- 
nalize T x and T so : Di = UmU^ 1 and D so = UsoT^U^ 1 . 
When one knows the matrix <f>i at an arbitrary site i, 
in principle, the wave function can be found anywhere 
by TM equations. We denote <F; = Ui$i in leads and 
$i = U so &i in middle bar. The TM equations in diago- 
nal representation can be written as 



A%, i G (-00, -L-1]U[L+ 1, 00), 
D so $i, 1 G [-L + 1.L-1]. 



(3) 



When the strip has no interface and fully homogenous 
or only one interface and two sides are homogeneous, it 
has been shown that the OTMM works well and has no 
numerical instability. In this letter, we study the strip 
with two interfaces. In this case, the OTMM will have 
serious overflow problem if the length of the middle bar 
between two leads is longer than a critical value. 

We denote D so and Di as D^ c \c = so, I. The di- 



k (c) 



D 



can be classified into 



agonal elements A,- 

two types: | = 1 and \\\ c ^\ ^ 1 which relate to 

propagating and evanescent modes respectively. For 



mode I A 



(c) 



1. it can be rewritten as e ±lkia , where 



fcj(> 0) is a real number, and a is the lattice constant, 
is called as "right going" wave, and e- lk ' a the "left 



a ikia 



going". For mode \\\°'\ 7^ 1, it can be rewritten as 
e ±('7i + l <W where r\i is a positive real number and fa 
the phase. The e r,ia+ ' l ^ i is a "right growing" or say "left 
decaying" mode, and e _,?ia_4 * i the "right decaying" 
or "left growing" . Due to det Ti = det = 1 in 
homogenous region, any modes e ±lkia or e ±W<*+*w 
must appear in pairs. Then we can always arrange 
the diagonal elements into the order: diag{D^ c '} 



{e 



ik±a 



p ik v a p — r}\a—Mp\ 



-iki a 



., e^ Q + 1 ^}. The first p + q states in diag{D} 
correspond to p right going and q right decaying modes, 
and second p + q states to p left going and q left decaying 
modes. At eigen energy E, totally we have 2p propagat- 
ing modes, and 2q evanescent modes. 2p + 2q = 4TV, TV 
is the number of chains. For each energy E, there are 
4TV modes distributed in 2N channels corresponding to 
2N different |Aj|, where |Ai.... p | = 1 and lAp+i^.^jvl > 1- 
If q equals zero, all states are extended. When p = 0, no 
propagating wave can exist in strip. Changing energy E 
results in changing of the p and q. 

We assume that a right going electron wave, e tklXi , 
is injected from the first channel of left lead. In 
general there should be some reflection waves in all 
channels in left lead, due to the scattering of inter- 
faces. We denote the reflection waves in 2N chan- 
nels as {r l e~ lklX \r p+m e r ' mXi+1 ' l,m : I = 1,2, ...,p;m = 
1,2, ...,q}. {rir p+m } describe the reflection coeffi- 
cients. The wave function in the left lead is $j — 
(e ifclXi , 0, 0, ne-* Ij , r p e- lk " x \ r p+1 e r ' lXi+i ' t ' 1 , 
r p+q e VqXi+ ' l '^ q ) t . We can further set the wave function at 



position i = —L to be $_l = (1, 0, 0, r±, r 2 , .., T2n) 1 ■ 
The phases of reflection waves have been absorbed in 
coefficients {ri}. We have $_l = A^-i-i, = 
T Sl $-l, $l = (^ ) 2L - 1 T Sl $-l, = 5$_ L ,where 



T s L - U so T s JJ l ,T Sr = UiTsjiU^ 1 

,2L-irp vv 



S = T Sr (D so ) 



Then we have = J2p=i Sa,p{$-L,)p = S a> i + 

Y^=i S a ,2N+ii"i) where a = 1, ...,47V. There are no 
left going waves injected from right lead, so the wave 
function at i = L + 1 can be expressed by $l+i = 
(ti, i 2 , t 2 N, 0, .., 0)*, where {ti} are the transmission 
coefficients. Thus we obtain 2N equations: 



(* i+1 ) a = 0,a = 2JV+l,- 



,47V. 



(5) 



Define 27V x 2N ^matrices W and Y : W im = 
S2N+i,2N+ m ;Yi = -S2N+i,i;l,m = 1,2, ...,2N and the 
reflection coefficient matrix R = (7*1, t 2 n) ■ Equa- 
tion ([5]) can be rewritten as WR = Y. The matrix R 
can be uniquely obtained through R — W _1 Y. Then 
27V transmission coefficients can be obtained from equa- 
tions t a = (4>L+i) a > a — 1,2,. ..,27V. Hence, the wave 
function at every site in leads and middle bar is ob- 
tained. So basically the electron transport though two 
interfaces and middle bar is completely solved. How- 
ever, the computation of W will encounter numer- 
ical overflow problem when L is large. The problem 
comes from the term (70 SO ) 2L_1 that may have some in- 
creasing modes with diagonal elements like |A p +i| -1 . 
If the value of |A p +i| 2i_1 3> 1, many elements Wi m — 

YA=i&S R )2N+Lt\i 2L ~ 1 (Ts L )i,2N+ m are of the order of 
lAp+jl Hence the calculation of W~ 1 will meet the 



numerical overflow with the order of {\ p+ i 2L ^ 1 ^j n , 1 <C 
n < 27V. Therefore, the OTMM can only accurately 
solve the solution for bar system with smaller length L. 
For example, we have calculated a Rashba SOC system 
with TV = 200 by means of OTMM. The longest length 
of bar is around 15, and numerical instability happens 
for 2L > 15. Larger TV results more shorter 2L. 

New improved TM method In fact, the physics here 
must be finite. The solution for {r^} should be stable. 
The numerical overflow is an artifact of the OTMM. Even 
\K\ 2L ~~ 1 ^S> 1 for evanescent modes when 2L — 1 ^> 1, 
the wave function after 2L — 1 steps of the transfer 
by (D so ) 2L 1 still should be finite. Thus the values of 
(f SL $- L ) 2 N+p+i,i = 1,2, . ..,<?, corresponding to the in- 
creasing modes of {|A p +i| > 1} must be very small to 

assure {{D so ) 2L 1 Ts L Q-L)2N+p+i to be finite. That is 
the physical requirement. Thus we introduce q new aux- 
iliary parameters {Q} and assume that 



(f SL ^ L ) 2N+p+i = Qe-^ L - 1 ^ ia+i ^,i=l,...,q. (6) 

q parameters {Q} have to be determined together with 
2TV reflection coefficients {r,}. The elements in matrix 
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<3>l+i is the linear combination of 2N coefficients {r^} 
and q auxiliary parameters {Q}- We have 



2 A 



^j=i 



(*£+l)a 

where 

= Z]0=l hP (^S , H)a/3(( £) so) 2i_1 ?S I ,)/31 



T.T=i P {Ts R U{{D so f L ^f SL )^ N+l , (8) 



— (Ts R )a:,2iV+p+i) a — 1,2, ...,4/V. 

No left going wave in right lead is requested, i.e. 
($L+i)a = 0,a = 2N + f, • • • ,4/V, which yield follow- 
ing 2iV equations 



L — *z= 1 * — '1=1 



(9) 



2iV equations ^ together with q auxiliary equations 
([t?|) , totally we have 2N + q equations that can uniquely 
determine 2N + q unknown parameters {r^} and {£?}• 
All coefficients of and {r^} are finite here so that the 
inversion of coefficient matrix can be calculated without 
overflow and unique solution of {Q} and {r^} are ob- 
tained. Now the numerical overflow problem does not 
exist any more. The method is stable for any length. Af- 
ter all {r"j, Q} are found, the transmission coefficients t a 
can be obtained by equations 



t a = ($ L+1 ) a ,a=l,2,...,2N. 



(10) 



As an example, we use our NITMM to study electron 
transport through a Rashba bar sandwiched by two semi- 
infinite metal leads. The Hamiltonian in bar region is 



H = 



Pi 



■n 



2m* 



a 



Px&y) +V coa {(x,y), (11) 



where (a x , a v , a z ) are Pauli matrices, a the strength of 
the Rashba SOC, and V con f(x,y) the transverse confin- 
ing potential. Here an open boundary condition in y 
direction is applied. TM in leads and Rashba bar can be 
written as following super-matrix 



/ A 

B* 



Ti 



B 
A 




V 



\ 




A 
B* 



B 

A) 



(12) 



where A and B are two 4x4 sub-matrix, A*(B*) is 
the complex conjugate matrix of A (£?). 



A = 



a b e f 

1 

— e — / a b 

10 



,B = 



g h 



-h* g* 





(13) 



where {a, 6, e, /, g, h} are the functions of the hopping 
constant to, SOC strength (t so in Rashba region, 
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(a) With EW 

E=-3.8t ,tso=0.02t ,200(w)x20(L) ' ™ 



(c) 
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FIG. 1: a), b), c) are the spin polarization(x=10) for 200x20, 
while d), e), f)(x=300) for 200x600, with the unit h/2. The 
continuous line included the evanescent waves(EW) . The dot- 
ted line included only the extending states. 



in lead region) and eigen energy E. The expres- 
sions of {a, 6, e, /, g, h} are a = —c{E — w x>y )/to,b = 
c (*so/*o - l) : e = -cbso (E - w x , y ) /*§, / = ~2ct so 
/t , g = c(it 2 so /tl - l), h = -c (1 - i) i so /t ,and c = 

(l + i 2 /to0 , where to = ft 2 /2m*a 2 , £ so = a/2a. tu^j, 
is the on-site energy, here we have chosen it to be zero. 

We take the E to be Fermi energy Ep, and fix the 
other parameters in TM, then consider an electron wave 
injecting from one of the channels in left lead and calcu- 
late the reflection and transmission coefficients {r^tj}. 
Finally the wave functions at every site of the strip can 
be obtained. We compare the result of our new improved 
method with that of the OTMM for 2L < 15 for N = 200 
system within which the calculation of OTMM is accu- 
rate and has no numerical instability. We obtain the 
exact same results for the transmission and reflection co- 
efficients. However, our new method can accurately cal- 
culate the case of 2L as large as we want. Main cost of 
computer time is from the diagonalization of 800 x 800 
TM. It is clearly shown that our new method does not 
require more computer time when the length of middle 
bar increases. Figure Q]and Ogive the exact solution of 
the spin polarization in middle bar and resonant trans- 
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FIG. 2: The upper orange line stands for transmission rate, 
while the lower black line is the reflection rate. 



algorithm in SMM. In ETMT0, they studied the mag- 
netoconductance of a quantum wire with several anti- 
dots. The numerical technique used there could avoid nu- 
merical instability happened in original transfer-matrix 
method, and was applied to a variety of 3D systems in- 
volving complicated atomic and many-body potentials. 
However, due to its iterative calculations, the computing 
time of ETMT also increases linearly with the increase in 
the length of quantum wire even for the homogeneous re- 
gion. Our method shows the superiority in treating long 
homogeneous transfer region for that it does not cost ex- 
tra computing time as the length of homogeneous scat- 
tering region becomes large. The computing time is in 
the zeroth order of homogenous region length L, O(L ). 
For the model we calculated here, our method is much 
faster than SMM for long length system, and we believe 
that it is also faster than ETMT. Our method can also 



mission of current as a function of bar length up to 10 
lattices, which is a size that can not be calculated by 
other methods. Where we have taken the average for 
non-polarized injection. The contribution from evanes- 
cent waves to spin polarization of Rashba SOC system 
based on continuous version has been studied by Li and 
Yangpjj]. From our results of figure [TJ the exact calcula- 
tion shows that the contribution from evanescent wave is 
significant for a short bar, but it is small when the length 
of bar becomes long. 

Comparison between our new method and previous 
methods (SMM and ETMT). SMM has been mostly ap- 
plied for studying the wave transport and has no problem 
of numerical instabilty [l3[ . In SMM, in order to avoid the 
numerical instability, one has to divide the middle bar 
into many small sections. Then one should find the SM 
for each section via TMM. In the sample with N = 100, 
the longest length of such pieces should be less than 400 
lattices such that no numerical instability happens. Here 
we use the 300 to ensure the precision of calculation. 
Then a recursion approach is applied to combine every 
sub SM to reach the final total SM[2l|. Fig.[3]shows that 
the CPU cost of SMM linearly increases with growing 
length of scattering region, but our new TMM is inde- 
pendent of the length of middle bar as showed in same 
figure. The numerical solutions obtained by SMM and 
our new TMM are exactly the same. In addition, when 
one studies the transmission of a polarized incident wave 
going through the middle bar, considering only a sin- 
gle channel injection of electron wave is not sufficient for 
SMM. All channels injection from left and right must be 
calculated individually to obtain the SM. Thus, the com- 
puting time to obtain the SM of first block near interface 
is almost 4N times of NITMM. Further more, computing 
time grows linearly for adding every block by recursive 
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FIG. 3: The time cost of SMM(orange upper) increases as 
the first order of lengths of scattering region, where as that 
of NITMM(black below) is the zeroth. The maximum length 
of a single block within which the SMM is applicable is less 
than 400a, for that there is a matrix inversion operation in 
the SMM formulae [2l[ . Here we choose the maximum length 
as 300a to ensure the precision of calculation. 

be applied to the studies on the transport under uniform 
magnetic field. Extensions of our NITMM to other prob- 
lems will be our future work. 
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